% This file analyzes the experiments intended to reveal the tensor monopole

%% process experiment data
r=2e6;
DQ_correction=[2,1]; % because we always divide by the full omega
start_avg=1; % if we start from the 1st or discard the 1st avg;

alpha_sweep=linspace(0,pi/2,9);
alpha_sweep=[alpha_sweep,alpha_sweep([5,7,5,4,5,5])]; % repeated measurements at a few alpha. In the final data processing these are averaged

%% start the analysis
flag = {[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,-1,0],[1,0,1],[1,0,-1],[0,1,1],[0,1,-1]};
flag_name = {'\Omega^\alpha','\Omega^\beta','\Omega^\phi','\Omega^{\alpha\beta}','\Omega^{\alpha\bar{\beta}}',...
    '\Omega^{\alpha\phi}','\Omega^{\alpha\bar{\phi}}','\Omega^{\beta\phi}','\Omega^{\beta\bar{\phi}}'};

beta=0;phi=0;
[rabi_estimate_save] = QGT_rabi_freq_detune(alpha_sweep,beta,phi,0,r,flag);

load('./data/FigS8.mat')

for hlp_alpha=6 % 6 corresponds to Fig S8-12
    alpha=alpha_sweep(hlp_alpha);
    for hlp_flg=1:length(flag)
        
        for hlp_omega=1:2 % omega   
            [rabi_fit,drabi_fit,pay_attention]=Rabi_for_tensor_monopole(x{hlp_omega,hlp_flg,hlp_alpha},y{hlp_omega,hlp_flg,hlp_alpha},dy{hlp_omega,hlp_flg,hlp_alpha},rabi_estimate_save(hlp_omega,hlp_flg,hlp_alpha)*E_pert_save(hlp_omega,hlp_flg,hlp_alpha),hlp_omega,flag_name{hlp_flg});
            
        end
%         saveas(gcf,['Rabi_data_',num2str(hlp_flg)],'epsc')
    end
end


function [rabi_fit,drabi_fit,pay_attention] = Rabi_for_tensor_monopole(x,y,dy,rabi_estimate,hlp_omega,name)
% a few principles: 1. have a relatively small range for contrast; 2.
% decide if t_period>>8us, so we do not count it's contribution; 3. start
% from rabi_estimate
% finally, check pbest vs delta_p. If on the same order, ditch the fit and
% give a warning

myfun_base =@(c,x) c(3)+ c(2)*cos(2*pi*x*c(1)+c(4));
amp=abs(x(floor(length(x)/2))-x(1));

switch hlp_omega
    case 1 % SQ
        myfun = @(c,x) myfun_base([c,pi],x);
        LB=[0.7*rabi_estimate,0,-0.3];
        UB=[2*rabi_estimate,0.1,0.05];
        pinit=[rabi_estimate,0.065,-0.215];
        
    case 2 % DQ
        myfun = @(c,x) myfun_base([c,0],x);
        LB=[0.7*rabi_estimate,0,-0.4];
        UB=[2*rabi_estimate,0.2,0.05];
        pinit=[rabi_estimate,0.13,-0.13];
end

if rabi_estimate<1
    pinit(end)=mean(y);
end

[pbest,delta_p]=easyfit(x, y, pinit, myfun, LB, UB,'np');

x_fit=linspace(0,9e-6,91);

if hlp_omega==1
    figure
end
cc=colormap(piratepal(9,1));
hold on
plot(x_fit*1e6,myfun(pbest,x_fit),'LineWidth',1.5,'Color',cc(hlp_omega,:),'HandleVisibility','off')
errorbar(x*1e6,y,dy,'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp_omega,:));

if hlp_omega==2
    xlabel('Duration of parametric modulation (\mus)','FontSize',18)
    ylabel('PL signal (a.u.)','FontSize',18)
    set(gca,'FontSize',18)
    lgd=legend({['SQ: $',name,'$'],['DQ: $',name,'$']},'FontSize',18);
    set(lgd,'Interpreter','latex');lgd.FontSize=18;
    ylim([-0.33,0.03])
end

rabi_fit=pbest(1);

if pbest(1)<1 % rabi frequency=0
    rabi_fit=0;
    drabi_fit=0;
else
    drabi_fit=delta_p(1);
end

if abs(rabi_fit-rabi_estimate)>0.05e6
    pay_attention=1;
else
    pay_attention=0;
end
end

function [rabi_g,E_total,full_map,g_munu_save] = QGT_rabi_freq_detune(alpha_sweep,beta,phi,rd,r,flag)

rabi_g = zeros(2,length(flag),length(alpha_sweep));
E_total = zeros(2,length(alpha_sweep)); % 1: SQ, 2: DQ
full_map = zeros(1,length(alpha_sweep));
dpert=0.005;
g_out = zeros(length(flag),length(alpha_sweep));
g_munu_save = zeros(6,length(alpha_sweep));

for hlp=1:length(alpha_sweep)
    alpha=alpha_sweep(hlp);
    H0=Ham_detune(alpha,beta,phi,rd);
    [V,E]=eig(H0);
    [E,idx]=sort(diag(E));
    
    E_total(:,hlp)=[abs(E(2)-E(1))*2;abs(E(3)-E(1))]*r; % note we already take 2X the SQ frequency here
    
    V=V(:,idx);
    um=V(:,1)/sqrt(V(:,1)'*V(:,1));
    u0=V(:,2)/sqrt(V(:,2)'*V(:,2));
    up=V(:,3)/sqrt(V(:,3)'*V(:,3));
    
    for hlp_flg=1:length(flag)
        flg_tmp=flag{hlp_flg};
        
        alpha_oscillate=flg_tmp(1);
        beta_oscillate=flg_tmp(2);
        phi_oscillate=flg_tmp(3);
        
        for hlp_omega=1:2 % omega
            dH=Ham_detune(alpha+dpert*alpha_oscillate,beta+dpert*beta_oscillate,phi+dpert*phi_oscillate,rd);
            dH=(dH-H0)/dpert;
            
            switch hlp_omega
                case 1
                    rabi_g(hlp_omega,hlp_flg,hlp)=abs(um'*dH*u0*r);
                    
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        abs(um'*dH*u0*r)^2/((E(1)-E(2))*r)^2;
                    
                case 2
                    rabi_g(hlp_omega,hlp_flg,hlp)=abs(um'*dH*up*r);
                    
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        (rabi_g(hlp_omega,hlp_flg,hlp))^2/((E(1)-E(3))*r)^2;
            end
        end
    end
    
    g_munu_save(1,hlp)=g_out(1,hlp);
    g_munu_save(2,hlp)=g_out(2,hlp);
    g_munu_save(3,hlp)=g_out(3,hlp);
    g_munu_save(4,hlp)=(g_out(4,hlp)-g_out(5,hlp))/4;
    g_munu_save(5,hlp)=(g_out(6,hlp)-g_out(7,hlp))/4;
    g_munu_save(6,hlp)=(g_out(8,hlp)-g_out(9,hlp))/4;
    
end
end

function Ham_out = Ham_detune(alpha,beta,phi,rd) %
Ham_out=[0,cos(alpha)*exp(-1i*beta),0;cos(alpha)*exp(1i*beta),0,sin(alpha)*exp(1i*phi);...
    0,sin(alpha)*exp(-1i*phi),0]+diag([rd,0,-rd]);
end